* Check correlation of OOH with all hours results - creates Figure A4

* Build estimates for OOH
use "$savedata/masterdata_09onwards.dta", clear

drop if exittime==.
keep if dow==0 | dow==6 | exittime<800 | exittime>2000

egen vol_trust = sum(one), by(pconsult trust_num)

foreach x in 10 25 50 100{
egen vol`x' = sum(one) if akm_set_mainspef`x'==1, by(pconsult)	
gen sample`x'=0
replace sample`x'=1 if akm_set_mainspef`x'==1 & vol_trust>=`x'
}

keep if sample25==1
gen vol = vol25

reghdfe survive30 c.prevyear_cost, absorb(i.sex##i.derv_age i.black i.mixed i.chinese i.asian i.race_miss i.ynch* i.prevyear_stroke i.di1 i.di2 i.di3 i.di4 i.di5 i.shock i.arythmia i.arthero i.arrest i.dow##i.admidate_mont##i.finyear i.hyid, savefe) keepsingleton
predict residual, residuals

xtreg residual c.std_ami3 c.std_nonami3, fe i(doctor_id)
predict docfe, u
gen sigma_u = e(sigma_u)
gen sigma_e = e(sigma_e)

collapse (mean) docfe residual sigma* vol, by(doctor_id)

gen signal_2step = ((sigma_u^2) / ((sigma_u^2) + ((sigma_e^2)/vol)))
gen adj_docfe = docfe*signal_2step

save "$savedata/09ooh_estimates.dta", replace

* Then produce main estimates
use "$savedata/masterdata_09onwards.dta", clear

drop if exittime==.

egen vol_trust = sum(one), by(pconsult trust_num)

foreach x in 10 25 50 100{
egen vol`x' = sum(one) if akm_set_mainspef`x'==1, by(pconsult)	
gen sample`x'=0
replace sample`x'=1 if akm_set_mainspef`x'==1 & vol_trust>=`x'
}

keep if sample25==1
gen vol = vol25

reghdfe survive30 c.prevyear_cost, absorb(i.sex##i.derv_age i.black i.mixed i.chinese i.asian i.race_miss i.ynch* i.prevyear_stroke i.di1 i.di2 i.di3 i.di4 i.di5 i.shock i.arythmia i.arthero i.arrest i.dow##i.admidate_mont##i.finyear i.hyid, savefe) keepsingleton
predict residual, residuals

xtreg residual c.std_ami3 c.std_nonami3, fe i(doctor_id)
predict docfe, u
gen sigma_u = e(sigma_u)
gen sigma_e = e(sigma_e)

collapse (mean) docfe residual sigma* vol, by(doctor_id)

gen signal_2step = ((sigma_u^2) / ((sigma_u^2) + ((sigma_e^2)/vol)))
gen adj_docfe = docfe*signal_2step

rename docfe docfe_all
rename adj_docfe adj_docfe_all

merge 1:1 doctor_id using "$savedata/09ooh_estimates.dta"

corr adj_docfe*
corr docfe*

twoway scatter adj_docfe_all adj_docfe, mcolor(gs3) || lfit adj_docfe_all adj_docfe, ytitle(Estimated doctor quality (All hours)) xtitle(Estimated doctor quality (Out of hours)) graphregion(color(white)) ///
legend(off)

graph export "$results/FigureA4.pdf", as(pdf) replace
